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ABSTRACT 

Performing a stable, long duration simulation of driven MHD turbulence with a high thermal Mach 
number and a strong initial magnetic field is a challenge to high-order Godunov ideal MHD schemes 
because of the difficulty in guaranteeing positivity of the density and pressure. We have implemented a 
robust combination of reconstruction schemes, Riemann solvers, limiters, and Constrained Transport 
EMF averaging schemes that can meet this challenge, and using this strategy, we have developed a new 
Adaptive Mesh Refinement (AMR) MHD module of the 0RI0N2 code. We investigate the effects of 
AMR on several statistical properties of a turbulent ideal MHD system with a thermal Mach number 
of 10 and a plasma f3o of 0.1 as initial conditions; our code is shown to be stable for simulations with 
higher Mach numbers (A^rms = 17.3) and smaller plasma beta (/3o = 0.0067) as well. Our results show 
that the quality of the turbulence simulation is generally related to the volume-averaged refinement. 
Our AMR simulations show that the turbulent dissipation coefficient for supersonic MHD turbulence 
is about 0.5, in agreement with unigrid simulations. 

Subject headings: Magnetic fields — MHD — ISM: magnetic fields — ISM: kinematics and dynamics — 
stars:formation — methods: numerical — turbulence 



1. INTRODUCTION 

Eulerian codes are commonly used in star formation 
studies in order to model the complex physical processes 
involved, including turbulence, magnetic fields, gravita- 
tional collapse, and radiation feedback. The dynamic 
ranges of density and size scales involved in star forma- 
tion are enormous, ranging from more than 10 pc in giant 
molecular clouds (GMCs) of dens ity ~ 5 AU protostella r 
cores with densities > 10^^ cm~^ (jMasunaga et al.iri998l ). 
This poses a significant challenge to numerical simula- 
tions using a uniform computation al mesh. For example, 
using the unigrid code ZEUS-MP ()Haves et al.ii2006l) . a 
simulation of ideal MHD turbulence with a 1024^ grid re- 
quires ~ 50,000 cpu hours. When gravitational collapse 
begins, dense cores will re ach the numerical resolution 
limit (Truelo ve et al.l |1997( ) in just a small fraction of 
the global free-fall time, tg. For a high-order Godunov 
scheme, the computing time will be ~ 5 times that for 
ZEUS-MP, which uses a low-order scheme. The comput- 
ing time is further increased by a factor of at least 16 
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whenever the resolution of the 3D grid is doubled be- 
cause maximum Alfven speed will be increased at lower 
density regions as the result of higher resolution. It is 
computationally inefficient to simply increase the grid 
resolution for star formation simulations because only a 
small fraction of the simulated region has collapsed to 
sufficiently high density to violate numerical resolution 
requirements; most of the simulated volume is in low 
density voids where such fine resolution is unnecessary. 
Therefore, adaptive mesh refinement (AMR) becomes an 
important tool for simulation of star formation using Eu- 
lerian codes. With AMR, the computational mesh is re- 
fined only in the localized regions where high resolution 
is required, and as a result computational resources are 
concentrated in the regions where they are needed most. 

Stars form in molecular clouds, which are cold (T ~ 
10 K), supersonically turbulent (sonic Mach numbers 
~ 10 on scales ~ 10 pc), and magnetized {B > lO^G); 
the Alfven Mach number is observed to be of o rder unity 
and the plasma (3 parameter is small (< 0.1; iCru tchei ] 
11999( 1. There are several MHD codes with AMR ca- 
pability, includ ing the pu blicly available codes Ramses 
(|Tevssieri[200l PLU TO (IMignone et al.l \2(M\ . ENZO 
(IWang fc Abef [20091) . and FLASH (IFrvxell et al1[2000[l . 
However, to our knowledge, there is no AMR code in the 
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literature that has demonstrated the capability of simu- 
lating supersonic MHD turbulence with initial conditions 
appropriate for star-forming regions. The primary rea- 
son for this is that Godunov schemes for ideal MHD with 
high-order approximate Riemann solvers cannot guar- 
antee posit hdW in density and pressure fe.g. |LeVeaud 
IT99I iTbrolfTOOa IBerthonll2005l: IZhang fc Shull2010n and 
are therefore unstable for turbulence that is driven for 
long times with such initial conditions. This becomes an 
important consideration when developing an AMR ideal 
MHD code for star formation simulations with MHD tur- 
bulence. 

Because turbulence is intermittent, one might hope 
that AMR would be very effective in simulating it. How- 
ever, in a study of p urely hydrodyn amic turbulence us- 
ing the ENZO code, iKritsuk et all (2006) argued that 
the use of AMR only became practical from an efficiency 
standpoint if the base mesh had high resolution to start 
with. When they attempted runs with coarse (128^ and 
256^) base meshes (note that they used refinement ratios 
of 4), a large fraction of their domain was refined until 
they were sufficiently able to resolve the localized turbu- 
lent structures using AMR: with the refinement criterion 
they adopted, they refined 90% of the computational do- 
main at 512'^ resolution, 65% of the domain at 1024'^ 
resolution and 34% at 2048'^ resolution. Furthermore, it 
is not just a matter of whether AMR is economical for 
a turbulence simulation, but also whether it can accu- 
rately capture the properties of the turbulence. Another 
attempt at using different refinement criteria, such as 
refining on local vorticity and divergence of velocity in 
addition to shock refinement for purely hydrodynamic 
turbulence, also shows that very large refinement cover- 
age generally result s in capturing the turbulence statis- 
tics HSchmidt et"aIll2009D . 

In this paper, we present a robust MHD AMR scheme 
that is able to simulate turbulent flows with high Mach 
numbers and strong initial magnetic flelds and that 
uses an accurate CT scheme to maintain V • B = 
to machine accuracy without recourse to approximate 
methods that rely on either divergence cleaning(e.g. 
ICrockett et al.|[2005l ) or monopole ad vection (Powell's 8- 
wave scheme and i t s extensions: e.g. IPowell et al.|[T999l: 
Dedner et al.' [ 20021: pyTignone et al.l I2007t iWang fc Abell 
2009; Waagan'sFalJ 120111). Recently, [Waagan (2009) 



modified the MHD module of the FLASH code using a di- 
rectionally split MUSCL-Hancock scheme with properly 
discretized Powell source terms in order to enable stable 
driven-t urbulence simu l ations . The driven-turbulence 
tests in I Waagan et al.l (j2011[ ) are either at high Mach 
numbers with a relatively weak initial magnetic field 
(/3o ^ 1) or at low thermal Mach numbers (< 2) with 
somewhat stronger fields (/3o ~ 0.25). The tests they 
carried out were all unigrid: the performance of their 
code on driven turbu lence with AMR was not described. 
iWaagan et al.l ()2011l ) cites simulations of driven turbu- 
lence at high Mach numbers in strong fields using this 
code, but these too were unigrid simulations. 

Here we present the results of a quantitative inves- 
tigation of simulations of ideal MHD turbulence using 
a newly implemented ideal MHD module in our AMR 
code, 0RI0N2, which is based on a conservative high- 
order Godunov scheme. We investigate the significance 



of refinement coverage to the quality of turbulence statis- 
tics in strongly supersonic MHD turbulent flow. We are 
able to perform long duration, high Mach number, driven 
MHD turbulence simulations with a magnetic field that 
is initially moderately strong {j3o ~ 0.1). In §2, we 
briefiy describe the numerical method and the implemen- 
tation of our AMR constrained transport scheme using 
the Chombo AMR framework ()Colella et al.ll2000[ ). In 
§3, we present several standard tests to examine the ac- 
curacy of the code for both unigrid and AMR simula- 
tions. In §4 we discuss the effects of refinement coverage 
on several ideal MHD turbulence statistics using AMR. 
We focus our discussion on the velocity power spectrum, 
the PDF, and the turbulence dissipation rates. In §5 we 
present our conclusions. 

2. NUMERICAL METHOD 

2.1. 0RI0N2 

In this paper we present 0RI0N2, which represents 
a major upgrade of ou r parallel radiation hydrodynamic 
AMR code "ORION" CTruel ove et "all[l99l iKlein IIT99I 
Krumholz et al. 2004, 2007). 0RI0 N2 is impleinented 
using the Chombo AMR framework (iColella et al.ll2000( ) 
and is in a modular form that allows the selection of 
a variety of physics modules for simulations. Chombo 
is a set of highly optimized tools that provide an in- 
frastructure for implementing finite difference and finite 
volume methods for solving partial differential equations 
on a block-structured AMR grid configuration. Elliptic 
and time-dependent modules are included, as well as sup- 
port for parallel platforms using MPI. The newly devel- 
oped MHD module of 0RI0N2 is based on the Godunov 
scheme in the finite volume formalism implemented in 
the PLUTO code (Migno nc e t al. 2007). Specifically, we 
use a dimension ally unsplit C orner Transport Upwind 
(CTU) scheme ([Colellal |1990D incorp orating the Con- 
strain ed Transport (CT) framework (jEvans fc HawlevI 
19881). Thi s CTU -HCT integrator is best described in 



Stone et al.l ()2008| ). which will not be repeated here. 



Our implementation preserves some of the fiexibility of 
the PLUTO code in choosing (1) different interpolation 
schemes to re-construct cell interface states from the 
cell-center values, such as the piecewise linear method 
(PLM) or the piecewise parabolic method (PPM); (2) 
different limiters for the preservation of monotonicity 
near a discontinuity during the re-construction stage, 
from the very diffusive minmod to least diffusive mono- 
tonized central difference limiter; (3) different Riemann 
solvers, such as Roe and HLL-family solvers, to obtain 
fluxes at the cell interfaces based on the reconstructed 
cell interface states; and (4) different CT electromotive 
force (EMF) averaging schemes, such as the simple arith- 
metic averaging of flux es computed during the upwind 
step (jBalsara fc Spiceil [19991. or a face-to-edge integra- 
tion procedure using the arithmetic average of the EMF 
derivatives from neighbor cells, or selecting EMF deriva- 
tives accordin g to the sign of the mas s flux at the cell 
interface. See [Gardiner fc Stond (|2005f ) on how to com- 
pute the EMF at ce ll edges during the upwind step. 
iMignone et al.l (|2007[ ) have summarized the advantages 
and disadvantages of the combinations of different solvers 
and integration schemes. 
The CT scheme increases the complexity of the algo- 
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rithm, especially for an AMR code, and it requires ad- 
ditional memory for storing the face-centered magnetic 
field. However, the benefit of the CT scheme is that the 
solenoidal constraint V • B = is ensured to machine ac- 
curacy. For cell-centered MHD algorithms , divergence- 
cleaning methods (e.g. Crockett et al.l [20051 ) are used to 
ensure the solenoidal constraint, but they cannot guar- 
antee positivity of pressure (or energy) and therefore 
reduce the robustness of the code. The iDedner et al.l 
(|2002l ) app roach, which is used in some ce ll-centered 
codes (e.g. iPeng k AbeTI 120091 : iMignone fc' Tzcfcraco^ 
[2OT0h . allows magnetic monopoles to decay with time 
as they are transported to the domain boundaries. Some 
studies have reported that monopoles could lead to in- 
correct ju mp conditions and oth er spurious dynamical 
effects fe.g lBalsara fc SpicedHool IToth 2000). The cell- 
centered field is readily calculated with the CT scheme: 
once the face-centered magnetic field is updated, the cell- 
centered magnetic field is computed from a simple aver- 
aging of the face-centered magnetic fields. 

2.2. AMR Implementation 

Extending any uniform-mesh algorithm to AMR re- 
quires several modifications, primarily involving the cou- 
pling between the solutions at different resolutions. We 
follow the block-structure d AMR approach outlined by 
Berger and Colellal (|1989( ) and extended to MHD by 
Balsaral ^OOW . 

We begin with a uniform mesh that spans the domain 
and is denoted as AMR level £ = 0, with mesh spacing 
h^. (Superscripts on the symbols riref, t and At indi- 
cate the level of refinement, not a power.) When refine- 
ment is triggered by some criterion, such as a steep den- 
sity, velocity or magnetic field gradient, refined grids are 
constructed in logically-rectangular patches with mesh 
spacing — h°/n'^^f, which are grouped into AMR level 
1; here the refinement ratio nj^^j is some power of 2. If 
further refinement is desired, more AMR levels are con- 
structed. Each AMR level £ has a uniform mesh spacing 
= h^~^ /n-^cf^ ■ In general, n^^^ can have a non- uniform 
dependence on £, but in this paper we adopt ?^f<^^ = 2. 
Patches with the same resolution are organized into lev- 
els, which are then organized into a hierarchy of AMR 
levels. 

Fbllowing the approach outlined bv lBerger and Colellal 
(|1989( ). we refine in time as well as space, commonly 
known as "subcycling" - the solution on each AMR level 
is updated using a timestep At^ = At^~^/n^^^. If we 
know the solution on the entire AMR hierarchy at time 
t°, then we begin the multilevel update of the solution 
on all levels by updating the base level solution by At°, 
without regard for any finer levels. We then recursively 
advance any refined levels until the solution on the en- 
tire AMR hierarchy has reached t° -I- Ai°. Whenever 
the solutions on different levels reach the same solution 
time, they are "synchronized" , to ensure that the com- 
posite solution remains conservative and maintains the 
solenoidal nature of the magnetic field. 

In general, the update of the solution on each logically- 
rectangular patch follows the same approach as that for a 
uniform-mesh implementation. AMR-specific implemen- 
tation details fall into the following categories: 

1) Boundary conditions: Following common prac- 



tice, boundary conditions at the edge of a rectan- 
gular patch are handled by adding a ring of ghost 
cells around the patch sufficient to complete the 
stencils used to update the solution on the inte- 
rior of the patch (otherwise known as the "valid 
region"). The algorithm in this work requires 4 
ghost cells. Solution values in ghost cells (and the 
magnetic fields on the associated "ghost faces" ) are 
filled depending on the type of boundary they are 
associated with: 

1.1) Physical domain boundaries - if the patch 
boundary abuts a (non-periodic) physical do- 
main boundary, ghost cell values are set using 
standard discretizations of the relevant physi- 
cal domain boundary condition (e.g. Dirichlet 
or Neumann boundary conditions). 

1.2) Copy boundaries - if the patch is adjacent to 
other patches on the same level, ghost-cell val- 
ues are filled by copying interior (valid-region) 
values from adjacent patches. 

1.3) Coarse-fine boundaries - if the ghost cells 
are along a coarse-fine boundary between two 
AMR levels, coarse-level solution values are 
interpolated to fill in ghost regions. If the 
fine level is not at the same solution time 
t as the coarse-level solution, we interpolate 
the coarse-level data linearly in time, since 
t^Td ^ ^ *lei- Cell-centered conserved 
quantities are filled using piecewise linear in- 
terpolation that is limited to prevent the cre- 
ation of new maxima an d minima, following 
iBerger and Colellal H1989D . The face-centered 
magnetic field is interpolated by first per- 
forming a piecewise-linear interpolation of the 
coarse-level field onto the fine-level faces that 
overlie the coarse-level faces, and then linearly 
interpolating these interpolated values in the 
face-normal direction to fill the fine-level faces 
that do not overlie coarse-level faces. Note 
that this differs from the diver gence- fre e inter - 
polation scheme presented bv iBalsaral ()2001[ ) 
in that interpolated values are not necessarily 
divergence free. We have not found it nec- 
essary to ensure that the magnetic fields in 
ghost regions are divergence-free, because the 
interpolated ghost values are used only in the 
reconstruction scheme to compute the fluxes 
in the divergence-preserving CT scheme and 
are never directly used to increment the mag- 
netic fields themselves. 

2) Interpolation to newly-refined regions: As 

the solution evolves, the refined regions evolve with 
it through a regridding process. Previously re- 
fined regions are de-refined when finer resolution 
is no longer needed, while coarse-level regions are 
refined when finer resolution is required. In the de- 
refinement case, we simply average the fine-level 
solution to the newly-exposed coarse-level mesh. 
Newly-refined regions are filled by interpolating 
the coarse-level solution. Cell-centered conserved 
quantities are interpolated using piecewise linear 
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interpolation that is limited to prevent the cre- 
ation of new solution maxima and minima. The 
face-centered magnetic field is defined using a two- 
step process. First, we interpolate the coarse-level 
B field onto the new faces using the face-centered 
interpolation used to fill ghost faces at coarse-fine 
interfaces. Then, the newly interpolated values are 
projected using a v ariant of the face- c enter ed pro- 
jection described in lMartin &: Colellal ()200ClD to en- 
sure that they are discretely divergence-free. 

3) Synchronization: When the solutions on two 
AMR levels reach the same time, a series of syn- 
chronization operations is performed. First, the 
fine-level solution is averaged onto the covered re- 
gions of the coarser level. This includes both cell- 
centered conserved variables and the face-centered 
B. The coarse- and fine- level solutions have been 
updated using fiuxes that have been computed in- 
dependently, likely resulting in a loss of conserva- 
tion at coarse-fine interfaces. Conservation is main- 
tained throug h a flux-correct i on ste p similar to that 
used by Bcrg er and Colellal (|1989t ). which ensures 
that the same fluxes are used to update the coarse- 
and fine-level solutions across coarse- fine interfaces. 

Similarly, the magnetic field is unlikely to be 
divergence-free at coarse-fine interfaces because 
coarse- and fine-level magnetic fields have been up- 
dated independently. Following iBalsaral (|2001l ). we 
ensure a divergence-free B at coarse-fine interfaces 
through a correction step that ensures that the 
same electric field values are used to update the 
coarse- and fine-level B field adjacent to coarse-fine 
interfaces. 

3. STANDARD TEST RESULTS 

We have performed many tests of the 0RI0N2 MHD 
module, including the well-k nown standard tests, such 
as the EM wave families (e.g. iCrockett et alJ | 2005f). the 
iRvu fc Jonei (\1995) shock-tube tests, the iBrio fc Wul 
( 19881) shock-tube test, the fiel d-loop advection test, an d 
the MHD blast- wave test as in lGardiner fc Stond ()2005[ ). 
We present the results of only two of the shock-tube tests 
and the field-loop advection test here to demonstrate the 
second-order accuracy of the MHD module. 

3.1. Shock-tube Tests 

iRvu fc Jones! (|1995l ) developed a suite of shock- 
tube tests that are commonly used for testing 
MHD algorithms. In Figure [l] we present the 
results of just one of the tests with the set- 
ting of initial conditions {p,Vx,Vy,Vz, B^, By, B^, P) = 
(l,10,0,0,5/(47r)i/^5/(47r)i/2, 0,20) on the left and 
(1, -10, 0, 0, 5/(47r)i/2, 5/(47r)i/2, Q, l) on the right of the 
contact discontinuity; here p is the density, v the veloc- 
ity and P the gas pressure. The contact discontinuity is 
located at the middle of the shock tube, and we set the 
adiabatic index of the gas to be 7 = 5/3. The length 
of the shock tube is 1 and is resolved by 512 cells along 
the a;-direction. The initial condition on the velocity is a 
colliding flow with a magnetic field at an angle of 45° in 
the X — y plane. The results at a time 0.08 in code units 
(corresponding to 0.46 sound crossing times at the left 



side of the shock tube initially) are shown in Figure [1] 
and they agree w ith the magnitudes a nd locations of the 
shocks found by IRvu fc Jones! (|1995f ) to almost within 
the thickness of the lines. 

In Figure [21 we show the results of an other com- 
monly used shock-tube problem, that of !Brio fc Wu! 
(|1988f ) . The initial conditions of this test are 
ip,v^,Vy,v,,B,,By,B,,P) = (1,0,0,0,0.75, 1,0,1) on 
the left and (0.125,0,0,0,0.75,-1,0,0.1) on the right; 
here P is the gas pressure. The contact discontinuity is 
located at the middle of the shock tube and 7 = 2. The 
length of the shock tube is 1 and is resolved by 800 cells 
along the x-direction. The gas has no movement initially. 
The j/-component of the magnetic field changes sign at 
the contact discontinuity and has a pressure jump of 10. 
Figure [2] shows the results at a time of 0.1 (corresponding 
to 0.14 sound crossing times at the left side of the shock 
tube initially). We can see the compound wave, which 
is composed of an Alfven wave and a slow wave, in the 
figure. 

3.2. Field-Loop Advection 

Since !Gardiner fc Stond (|2005| ) first suggested using 
the advection of a magnetic field loop to show the differ- 
ence between operator split and unsplit schemes, field- 
loop advection has become a standard test for MHD al- 
gorit hms. Our setup of the 3D test is exactly the same 
as in [Gardiner fc Stoni ()2008D : the magnetic field loop 
is created inside a rectangular box of size (2,1,1) resolved 
on a 2A^ xNxN grid with periodic boundaries. The loop 
of magnetic field is generated from a vector potential on 
a coordinate system (xi,X2,X3), with Ai — A2 — and 

, _(BoiR~r) for r<i?, . . 

^3 - I for r > i?, 

where Bq — 10^^, r — y/xj+x^, and the size of the 
loop is -R = 0.3. The computational coordinate system 
{x,y,z) is transformed to {xi,X2,X3) by 

xi^ {-2x + z)/V5, 

X2=y, (2) 

X3~ {x + 2z)/V5, 

corresponding to a rotation about the y-axis. The den- 
sity (p = 1) and pressure (P — 1) are uniform and the 
whole region is advected with a velocity (2,1,1). There- 
fore, after one unit of time, a loop starting in the middle 
of the rectangular region will travel across the 3D diag- 
onal of the box and return back to the initial position. 
The advection continues for 2 cycles and the images of 
the magnetic field loop at the beginning and the end of 
this test are shown in Figure |3l The top part of the 
figure shows the simpler case in which the field loop is 
aligned with the z-axis, which is similar to a 2D advec- 
tion test. The loops diffuse slightly but maintain their 
shapes nicely after 2 cycles of advection. The time evo- 
lution of the volume mean SB^/B^ is similar to that 
of the 3D inclined field-loop test, and B^ remains zero 
to the machine accuracy. In the rest of this section, 
we focus on the 3D inclined field-loop test results. We 
have performed the inclined field-loop test 3 times on 
a single level grid (unigrid) with 3 different resolutions: 
TV = 32, 64 and 128, as in .Gardiner fc Stone. (,2008,) ■ No 
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AMR is used. The time evolution of the mean square 
field, {B"^), normalized by the initial value, {Bj), of the 
whole advection sequence is shown in Figure as the 
three thin curves. The results are similar to those of 
IGardiner fc Stonj (|2008| ). who used a second-orde r PLM 
recons truction scheme, and to those of iFromang et al.l 
(|2006[ ), who used a second-order TVD scheme for the 
Ramses code. IGardiner fc Stonj ()2008[ ) pointed out that 
with the axis of the field loop aligned along an inclined 
direction with respect to the computational grid, pre- 
serving i?3 to be zero is non-trivial. The time evolution 
of the normalized error, < \B^\ > /Bo, of the 3 tests is 
sh own in Figure Bb; the erro r is slightly smaller than that 
in IGardiner fc Stonj (|2008D . Figure El shows the volume 
rendering of the inclined field loop from the A'^ = 64 un- 
igrid model after 2-cycles of advection along the 3D box 
diagonal. 

We next investigate how AMR affects this test. In Fig- 
ure ini we show an AMR advection test with a base grid 
of resolution A^ = 64 and one level of refinement, with a 
refinement ratio n^cf = 2. The refinement criterion is on 
the jump in the normalized magnetic pressure {SB^ / B^), 
and the refinement threshold is 2.8. With this refinement 
criterion, the entire field loop is covered by the fine grid. 
The level 1 fine grids continually move with the loop. Af- 
ter 2 cycles, the loop maintains its shape as in the advec- 
tion test on a unigrid with a resolution of A^ = 128. The 
dotted curve in Figure |3^ almost exactly coincides with 
the thin solid curve. The normalized error {\B^\)/Bq is 
also close to that of A^ = 128 test. 

In Figure [71 we show an advection test using fixed mesh 
refinement (FMR) as opposed to adaptive mesh refine- 
ment (AMR). The level 1 fixed fine grid, which also has 
rii-of = 2, is smaller than the base grid in space, and there- 
fore the loop passes through the coarse-fine grid bound- 
ary during the advection. The loop maintains its shape 
nicely after 2 cycles of advection passing in and out of 
the fine grid, but it diffuses slightly after the first pas- 
sage through the coarse/fine boundary. The normalized 
values of (i?^) and (IB3I) are shown in Figure and 
lib (thick dashed curve). Note that this test also high- 
lights the need for effective refinement criteria and for 
adaptive refinement that follows the solution, since this 
test only results in comparable accuracy to the uniform 
coarse mesh. In the unigrid, AMR and FMR tests, V • B 
vanishes to order 10~^^ (machine accuracy). 

4. SUPERSONIC ISOTHERMAL MHD TURBULENCE WITH 

AMR 

4.1. Simulation Model Parameters 

In this section, we investigate the effects of AMR on 
ideal MHD simulations of isothermal, supersonic turbu- 
lence in a strong magnetic field. We discuss only the 
effects of AMR on the velocity power spectrum, the den- 
sity PDF, and the turbulent dissipation rate. We present 
the results of 10 simulations, all with an rms 3D sonic 
Mach number A^rms = 10 and an initial value of the 
plasma /3 paramater /3o = 0.1; the corresponding initiial 
Alfven Mach number is Ma,o = -v/S- Periodic boundary 
conditions are used. We implemented an algorithm for 
driving t he turbule n ce wi th AMR using the recipe dis- 
cussed in iMac Low! (|1999D . The high value of the sonic 
Mach number and the low value of /3o have proven to be 
very challenging in the past since the simulation can be- 



come unstable within or soon after one dynamical cross- 
ing time, especially with AMR. We define the dynamical 
crossing time as the length of the box divided by the rms 
velocity of the turbulent box. The driving pattern we 
use is the same for all 10 models to facilitate direct com- 
parison. The system is continuously driven at the largest 
scales, fc = 1 ^ 2, so as to maintain A^rms = 10 for 3 
dynamical crossing times (i.e., 0.3 sound crossing times). 
We analyze the data only after the first crossing time so 
as to allow the system to relax to a steady state. There 
are a total of 50 data files written out from each model 
in the last two dynamical crossing times. 

We have explored different interpolation schemes and 
Riemann solvers in order to determine a combination 
of these algorithms which provides accuracy and stabil- 
ity for driven turbulence over several dynamical crossing 
times. Here we describe the combination of algorithms 
we adopted for our tests: 

1) The second-order accurate in space, piecewise to- 
tal variation diminishin g (TVP) lin ear interpola- 
tion scheme described irT lTorol ()1999[ ). 

2) The multi-d imensional shoc k flattening strategy 
developed bv lMignond (j2005D . in which the interpo- 
lation reverts to the minmod limiter and the fluxes 
are computed using the HLL solver when a strong 
shock is detected. This provides additional dissi- 
pation in the proximity of a strong shock so as to 
guarantee positivity of the pressure. 

3) The harmonic mean limiter of Ivan Leeil ()1974D . 

4) With the CT scheme, the electromotive force 
(EMF) is computed at the zone edges using a two 
dimensional Rieman n solver based on a four-state 
HLL flux function (jlondrillo k del Zannal 120041: 
iMignone et al.ll2007l) . 

5) The simple three-state HLLD approximate Rie- 
mann so lv er for the isothermal case described by 
IMignone! ()2005D . The absence of the entropy 
mode in the isothermal case leads to a differ- 
ent formulation based on a three-state represen- 
tation rather than t he fou r-state representation of 
IMivoshi fc Kusanol (j2005D . The MHD module can 
also handle a non-isothermal ideal gas, but we do 
not include tests with non-isothermal turbulence 
here. 

6) The ch aracter istic tracing scheme of 

IColella fc WoodwaTdl (fT98l . 

We have tried the Roe solver (jRod I1986[) , which is 
an approximate linear Riemann solver, with the above 
combination, and it is equally stable for long-duration 
driven turbule n ce sim ulations. The standard tests in 
IMignone et~all ()2007[ ) show that the Roe and HLLD 
solvers yield comparable accuracy, but the HLLD solver 
is faster. Our tests lead to the same conclusion, so we 
use the HLLD solver for all the driven MHD turbulence 
simulations in this investigation. 

To push this stable scheme further, we carried out two 
simulations on a 128"^ base grid with 2 levels of refinement 
for 3 dynamical crossing times: (1) A^rms = 10, plasma 
Po = 0.02; and (2) TW^ms = 17.32, plasma ^0 = 0.00667. 



6 



Li, McKee, & Klein 



We found that the former test is stable using a CFL 
number of 0.4 and the latter test is stable using a CFL 
number of 0.35. The results for the second test (Model 
11) are shown with the other 10 models for reference, 
even though the initial conditions of this test are dif- 
ferent. All the standard tests shown in S|3]use the above 
combination of algorithms and demonstrate the accuracy 
resulting from this choice. However, it must be borne 
in mind that our choice of Riemann solver, reconstruc- 
tion, limiter, and EMF averaging schemes is probably not 
unique, since we have tested only a small combination of 
all existing solvers and reconstruction schemes. 

We have also tried the thi rd-order accurate Piecewise - 
Parabohc-Method (PPM) of lColella fc WoodwardI ([l98l 
as the reconstruction scheme, but the simulation be- 
comes unstable within or soon after one dynamical cross- 
ing time, regardless of what limiter is used (even with the 
most diffusive minmod limiter). We conclude that spe- 
cial treatments will be required to use higher-order inter- 
polation schemes for simulations of high Mach number, 
strong magnetic field turbulence. 

4.2. Refinement Criteria 

iKritsuk et al.l (|2006[ ) have proposed two refinement cri- 
teria for hydrodynamic turbulence, one based on the 
jump in pressure and one on the norm of the veloc- 
ity gradient matrix. We have used this as a guide- 
line in setting our refinement criteria for MHD turbu- 
lence. For the first criterion, we replace the thermal 
pressure P by the sum of the thermal and magnetic pres- 
sures, Ptot = -P + B^/8tt; cells are tagged for refinement 
if APtot /Ptot exceeds a pre-de termined threshold. For 
the second criterion we use the IKritsuk et al.l (|2006D re- 
finement criterion for strong shear, which is determined 
by computing the norm of the velocity gradient matrix 
without the contribution from the diagonal ele- 
ments. The norm is then normalized by c^/Ax and is 
tagged for refinement at a threshol d that is the same as 
that for the total pressure jump. IKritsuk et ah* (2006^) 
found that AMR results agreed well with unigrid results 
at the maximum resolution of the AMR for a pressure 
jump threshold of 2. With the inclusion of the shear ve- 
locity refinement criterion, they found that the pressure 
jump threshold could be raised to 3. 

We have carried out a series of tests to investigate 
the sensitivity to the refinement threshold. Unlike 
IKritsuk et a"Ll(|200l ). we find that a refinement threshold 
APtot /Ptot = 2 causes the base level to be completely re- 
fined in all our AMR models, most likely bec a use o f the 
inclusion of magnetic pressure. IKritsuk et al.l ()2006D use 
an AMR refinement ratio nf^^ = 4. For example, their 
1024^ AMR run has a base grid of 256^ with one level 
of refinement. All of our tests use nf^^ = 2, except for 
one with n^^f = 4 for comparison. With a refinement 
ratio of 2, two levels of refinement are needed to achieve 
a maximum resolution of 1024^ with a 256'^ base grid. 
The fraction of the volume covered by the first level of 
refinement (i.e., the "coverage") is independent of the 
refinement ratio. Thus, the coverage at level 1 with a 
refinement ratio of 2 (maximum refinement equivalent to 
512^) is the same as that at level 1 with a refinement ra- 
tio of 4 (1024^), but the number of cells is 8 times larger 
in the latter case. The coverage at level 2 with a refine- 



ment ratio of 2 (1024^) is smaller than that of level 1 
with a refinement ratio of 4, and the number of cells is 
correspondingly smaller. 

A simple way of characterizing the refinement of a mul- 
tilevel AMR calculation is the volume-averaged resolu- 
tion, 

{R)=J2R^y^■ (3) 

i 

where Ri is the ID resolution of level i and Vi is the 
fractional volume coverage for level i excluding higher 
levels. For example, for Model 4, the volume coverage of 
levels 0, 1 and 2 is = 0.35, 0.48, 0.17, respectively (the 
total level 1 coverage is 0.65 but Vi excludes the volume 
refined at level 2). For this model, (P) = 128 x 0.35 + 
256 X 0.48 + 512 x 0.17 = 255, which is very close to that 
of a 256'^ unigrid model. 

We have performed a large number of ideal MHD tur- 
bulent box experiments using different base grid sizes 
(128^ and 256^), different refinement criteria, different 
minimum block sizes, and different refinement ratios. We 
also have two unigrid turbulence models, at 128^^ and 
512^, for comparison. The unigrid 512^ model will serve 
as the reference for all AMR models presented in this 
paper. 

4.3. Power Spectrum 

In this section, we study how changes in resolution 
and refinement criteria affect the velocity power spec- 
trum, Pv{k) cx in terms of the power index, n, 
and the extent of the inertial range, fcmax- We fit the 
power spectra obtained from the 50 data dumps between 
1 and 3 dynamical crossing times between = 4 (to avoid 
the effects of driving) up to a value of k that increases 
by unity at each iteration. All the fitting results and 
the plots shown in section [J] are time-averaged results 
from the 50 data dumps over these two dynamical cross- 
ing times. (Since the turbulence remain s corr elated for 
about one dynamical crossing time fe.g. ILi e~ al. 2008), 
we assign an error to the mean value of P(fc) from the 
50 data dumps equal to the standard deviation divided 
by -^3.) As more points are added, the uncertainty in 
the slope decreases, and correspondingly so does the re- 
duced of the fit. However, the reduced begins to 
increase when the power spectrum turns over due to nu- 
merical dissipation; we define the point at which is a 
minimum as the upper end of the inertial range, Zcmax- 
In order to overcome the possibility that noise in the 
data could artificially lower fcmax, we omit the point at 
fcmax + 1 that led to the increase in and evaluate the 
reduced of the fit including the point at fcmax + 2; 
if the value of the reduced is less than the previous 
minimum, we set fcmax = fcmax (previous) -I- 2 and proceed 
with the iteration. We allow for the possibility that the 
noise fluctuation could be up to three cells wide. This 
procedure allows our estimate of the inertial range to 
extend beyond the bumps at fc ~ 9 that are apparent 
in the spectra in Figure [5] and are due to an artifact in 
the driving pattern. This method is conservative, but 
it can eliminate the impact of the bottleneck effect on 
determinin g the spectral index or the size of the inertial 
range fe.g. IVermal 120071 : IKritsuk et"all I2007D when the 
inertial range available for turbulent energy transfer is 
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small. However, we note that our MHD turbulence tests 
do not suffer from the bottleneck effect. 

In Table 1, we present the fitted results of the 10 ideal 
MHD turbulence models. The table includes the refine- 
ment coverage at each fine level. For example, in Model 
4, level £ = 1 has a volumetric coverage of 65% of the 
computational domain and level 2 has a volumetric cover- 
age of 17% of the computational domain (or 17/65 ~ 26% 
of level 1). In Figure [U the compensated power spectra 
of models 1, 3, 4, 6, and 10 are shown. We summarize 
the results in Figure [8] and Table 1 as follows: 

1) Models 1 and 10 are unigrid simulations and pro- 
vide a basis for evaluating the AMR models. The 
spectral index of Model 10 is n = 1.42 ± 0.02, con- 
sistent with other strong magnetic field supersonic 
turbulence simulations, which show that the power 
spectrurn is close to the Iroshnikov-Kraichnan spec- 
trum (jlroshnikovi 1 1 9631 : iKraichnanI 1 1 9651 . The low 
resolution Model 1 has a steeper spectral index, 
n = 1.75 ± 0.06. We repeated Model 10 (uni- 
grid 512'^ with an HLLD solver) using the Roe 
solver and obtained a power spectrum that agrees 
to within the uncertainty of fitting. 

2) Models 2 to 4 test the effect of changing the refine- 
ment threshold for the total pressure jump from 
3.25 to 2.5. As the threshold drops, the refinement 
coverage increases, and the spectral index slowly 
approaches that of Model 10. Correspondingly, 
the inertial ranges are significantly longer (2 times) 
than in Model 1 and appear to converge to that of 
Model 10 as the refinement coverage increases. 

3) Model 5 tests the effect of shear flow refinement 
on the AMR calculation. There is no noticeable 
increase in the accuracy of Model 5 compared to 
Model 4, presumably because the additional cri- 
terion did not significantly increase the average 
refinement. Model 8 has shear flow refinement, 
whereas Model 7 does not; however. Model 8 has 2 
levels of refinement compared to 1 for Model 7, so 
no inference on the effect of the shear flow refine- 
ment can be drawn. 

4) Comparison of models 5 and 6, and of models 4 
and 7, addresses the effects of refinement coverage 
on the overall improvement of the power spectrum. 
Model 6 has only 1 level of refinement but uses a 
refinement ratio of 4, equivalent to full coverage 
of level 1 at the resolution of level 2 in Model 5, 
which is a 2-level model using a refinement ratio of 
2. The average refinement in Model 6 is about 60% 
greater than in Model 4; correspondingly, there is 
a substantial improvement to the power spectrum 
and a modest increase in the size of the inertial 
range. However, the computational time for Model 
6 is about 3 times that of Model 5, a big price to 
pay for the improvement. The reason for the large 
increase in computing time is that each cell that 
is refined only to level 1 in the two-level run has 
8 times as many refined cells in the one-level run. 
An additional test of the effects of refinement cov- 
erage is provided by comparing Model 7 to Model 
4; the former uses a base grid of 256'^ and only 1 



level of refinement. Model 7 has an average refine- 
ment about 20% greater than Model 4. For this, 
one gets a significant improvement in the spectral 
index, but only a small increase in the size of the 
inertial range. The computation time for Model 7 
is similar to that for Model 6. 

5) Like Models 2-4, Models 8 and 9 address the im- 
plications of increased refinement coverage for im- 
provements in the power spectrum. Note that the 
level 2 resolution in these models is the same as 
that in a 1024^^ resolution simulation. Model 9 
has only a slightly lower threshold for the pressure 
jump than Model 8 (2.3 vs. 2.5), but this leads 
to almost a twofold increase in the average refine- 
ment. This increase in average refinement yields 
an increase in the inertial range, but has no signifi- 
cant effect on the spectral index. Even though the 
average refinement of Model 9 exceeds that in the 
unigrid Model 10, the latter appears to have the 
most accurate spectral index and the largest iner- 
tial range. Our finding that a unigrid simulation at 
512'^ resolution is superior to an AMR simulation 
with a maximum refineme nt of 1024 '^ is c onsistent 
with the conclusion of Kri tsuk et al.l ()2006il that a 
large base grid is required to obtain the advantages 
of AMR for simulations of turbulence. 

6) Model 11 has different initial conditions from the 
above 10 models. This model has A^rms = 17.32 
and plasma /3o = 0.00667; the corresponding 
Alfven Mach number is Ma,o — 1- This model is in 
many ways similar to Model 3 in Table 1 in terms 
of the power spectral index, the length of inertial 
range, and refinement coverage, although Model 3 
has very different initial conditions (A^rms — 10, 
(3o = 0.1, and Ma^o = V5). 

From the above summary, we can see that changes in re- 
finement criteria directly affect the average refinement, 
which in turn directly affects the quality of the turbu- 
lence power spectrum. 

In Figure ini we show the power spectra of the magnetic 
field for models 1, 3, 4, 6, and 10. The convergence 
behavior of the magnetic field power spectrum is similar 
to that of the velocity power spectra in Figure |S1 

A recent study of the effects of purel y solenoidal and 
purel y compressive turbulent driving (iF ederrath et al.] 
120101 ) shows that there are significant differences in the 
turbulence statistics between these two extreme driv- 
ing models. For example, the dispersion in the den- 
sity PDF with purely compressive driving can be 3 
times that resulting from purely solenoidal driving. For 
purely solenoidal driving, the solenoidal component of 
the power spectrum will dominate the dilatational com- 
ponent. This is reversed for purely compressive driv- 
ing. Reality probably is somewhere in between these 
two extreme cases. Our driving is purely solenoidal. We 
decompose the solenoidal (V • = 0) and dilatational 
(V X Vc = 0) components from the velocity, v = + Vc, 

by 

v,(k) = [k.v(k)]k (4) 
v,(k) = [kx v(k)] xk (5) 
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(|Lemaster fc Stonil2009D . We define the fraction of the 
dilatational component in the velocity power spectrum 



Y ik) = ^ 



(6) 



(jKritsuk et alj l2009al) . Figure [TO] shows the time- 
averaged velocity power spectrum of Model 10 with 
the solenoidal and dilatational components. The time- 
averaged Xc(fc) = 0.27±0.01, 0.26 ±0.01, and 0.26 ±0.02 
for Models 1, 4, and 10, respectively. It appears that 
Xc{k) is not sensitive to the resolution. For Model 
11, which has a smaller Alfven Mach number, Xc{k) = 
0.27 ± 0.02 is the same. The value of Xc{k) in an ideal 
MHD turbulence mod el with purely solenoidal forcing 
(IKritsuk et al.l (l2q09aD) is « 1/4, similar to our values. 
ILemaster fc Stong(|2009D measure the dilatational com- 
ponent of the kinetic energy (i.e., the density- weighted 
velocity power spectrum) and find that it is even smaller 
compared to the solenoidal component. 

4.4. Density PDF 

In Figure lllb . we show the probability density func- 
tions (PDFs) of density for models 1, 4, and 10 to 
demonstrate the effect of refinement on the density PDF 
in AMR simulations. The density PDF of the unigrid 
Model 1 at 128'^ resolution is narrower than the density 
PDF of the unigrid Model 10 at 512'^ resolution. A low 
resolution unigrid turbulence simulation cannot reach the 
highest and lowest densities attainable in a high resolu- 
tion unigrid simulation. As a result, the maximum wave 
speed in a low resolution simulation can be lower than in 
a high resolution simulation, since the minimum density 
is higher and the maximum Alfven velocity is most likely 
smaller. On the other hand, since the maximum density 
in a low resolution simulation is smaller, the clump mass 
function will be reduced at high densities; this could be 
problematic in simulations of star formation. 

AMR offers the best of both worlds: As shown in Fig- 
ure [11^, at low densities the AMR simulation (Model 4) 
is close to Model 1 and therefore does not have the very 
high Alfven velocities that appear in the high resolution 
Model 10. In simulations of star formation, the loss of 
resolution at low densities is not important. However, 
Model 4 has the same maximum resolution as Model 10, 
and the AMR enables it to track accurately the high- 
density portion of the density PDF, which is critical in 
simulations of star formation. 

In Figure [TTb . the density PDFs of models 2 and 3 
are plotted on top of the density PDF of Model 10. The 
high-density part of the PDF of Model 2 deviates more 
from that of Model 10 than that of Model 3 because of 
the very low level 2 coverage (1.3%) in Model 2. Model 3 
suggests that in order to have a good match to the high 
density part of the Model 10 PDF, the level 2 coverage 
must be > 10% for this problem. With the Roe solver, 
the PDF of a unigrid simulation extends to slightly lower 
densities than that for the same unigrid simulation using 
the HLLD solver, but the high-density parts of the PDF 
are almost the same. Therefore, not only is the Roe 
solver slower than HLLD solver per time step, but the 
time step in an MHD turbulence simulation will also be 
smaller. 



4.5. Turbulent Energy Dissipation Rate 

Both hydrodynamic and MHD turbulence simulations 
show that turbulence decays on the order of a dynami- 
cal crossing time. The dissipation rate of the turbulent 
energy is of order pv'^^^/L. Specifically, we write 



E = e 



P'"Ls 



(7) 



where Wrms is the density-weighted rms velocity of the 
turbulence and the integral length scale (fBatchclor 19531) 
for a compressible fiuid with a magnetic field is defined 

by 



Lint — 



k-^Etot{k)dk, 



(8) 



where E^otik) is the total energy power spectrum, in- 
cluding both kinetic energy and magi ietic energy. Su - 
personic MHD tu rbulen ce simulations (Mac Low (199§), 
iStone et al.l (fT998i) . and ILemaster fc Stone (.2009)) sug- 
gest that th e pro portionality constant e ^ 0.5. 
IKaneda et al.l ()2003D summarized the results of many 
pure hydrodynamic incompressible simulations and 
showed that, with their highest resolution (4096"^) simu- 
lation, e converged to about 0.4. 

The turbulent dissipation rate coefficient e for the 10 
models presented here is plotted as a function of the av- 
erage resolution in Figure[T^ The horizontal uncertainty 
bar is obtained from the standard deviation of the vari- 
ation of the refinement coverage during the simulation, 
usually a few percent of the refinement coverage. Mod- 
els 1 and 10 do not have horizontal error bars since they 
are unigrid models. Figure [T^] shows that the dissipa- 
tion rate is within the uncertainties of the 512^ unigrid 
simulation for (R) > 350. For the three AMR mod- 
els with (R) > 350, the mean e = 0.48 ± 0.02, which 
agrees well with the result from the 512'^ unigrid model, 
e = 0.48 ± 0.01, and is similar to the results from other 
unig rid simulations of ideal MHD turbulence, e ~ 0.5 
(e.g. ILemaster fc Stonell2009l ). In Figure IT^ we also plot 
the value of e for Model 11 (solid circle) for reference. 
Although this model has the same spectral index for the 
velocity power spectrum as Model 3, the value of e of this 
model is smaller. 

Models 6, 8, and 9 have similar turbulent dissipation 
rates as Model 10, and one might think that using a lower 
resolution base grid with AMR would have the advantage 
of reaching the dissipation rate obtained from a unigrid 
model with higher resolution. However, Table 1 shows 
that the CPU usage for these 3 models actually is higher 
than for the unigrid Model 10. This is another indication 
of the inefficiency of AMR for turbulence studies until a 
much larger base grid is used. 

4.6. Turbulent Magnetic Field 

For a magnetized turbulent system, the magnetic field 
strength will be enhanced as the result of field-line 
stretching. The maximum field-line stretching will be 
limited by the grid resolution, which determines the nu- 
merical diffusion of the magnetic field that we want to 
minimize. We plot the time-averaged change in mag- 
netic field strength, (|Ai3|), normalized by the initial 
mean field strength, Bq, versus the volume-averaged 
resolution (i?) in Figure [131 Since the initial Alfven 
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Mach number is modest (A^a.o = V^) for Models 1- 
10, the enhancement of the field strength is also modest, 
{\AB\) ^ (1.54±0.09)Bo from the mean value from Mod- 
els 6 to 10, which have overlapping error bars and which 
have (R) > 300. The value of {\AB\)/Bo = 0.53 for 
Model 11 is smaller since A4a,o = 1 is smaller. Only sys- 
tems with initially high Alfven Mach numbers can have 
a turbulent magnetic field many times greater than the 
mean field, and simulations of such systems require much 
higher resolut ion to have conv erged magnetic field statis- 
tics fe.g. Krit suk et al.ll2009al) . 

5. CONCLUSIONS AND DISCUSSIONS 

We have developed a robust MHD module for our new 
AMR code, 0RI0N2, and have demonstrated its ability 
to carry out accurate long-duration AMR simulations of 
highly supersonic turbulent fiows with strong magnetic 
fields (/3 <C 1, Ma ~ !)• Although unigrid simulations 
of such fiows have been published (e.g., Krits uk et al.l 
r2009hl , to our knowledge the AMR simulations of these 
flows presented herein are the first to appear in the litera- 
ture. Since observations suggest that GMCs have highly 
supersoni c flows with re latively strong magnetic fields 
([McKec fc Ostrikeri[2007l) and since AMR is essential for 
following gravitational collapse, this represents an impor- 
tant advance in our ability to study the conditions that 
lead to star formation. 0RI0N2 is able to do this be- 
cause the code is sufficiently flexible that one can easily 
experiment with different reconstruction schemes, Rie- 
mann solvers, CT EMF averaging schemes, and limiters 
to find the suitable combination that we have described. 

We have tested the accuracy of our c ode with sev- 
eral standard MHD tests, in cluding the IRvu fc JonesI 
(|1995[) shock tube test, the IBrio fc Wul ([19880 shock 
tube test, and the 3D current-loop test on unigrid, 
FMR, and AMR. The results we have presented here 
demonstrate that the CTU-I-CT algorithm performs ac- 
curately and works properly within the Chombo AMR 
framework. We found that the piecewise linear, spa- 
tially second-order, TVD scheme combined with a multi- 
dimensional sh ock flattening strateg y developed by the 
PLUTO group (jMignone et al.ll2007n can work with both 



Roe and HLLD solvers to enable stable, long-duration (3 
crossing times) simulations of driven MHD turbulence 
with an rms Mach number of A^ims = 17.3 and an initial 
plasma (3 parameter of /3o = 0.0067 on a 128"^ base grid 
with 2 levels of refinement. 

We examined the velocity power spectrum, the density 
PDF, and turbulent dissipation rate in our investigation 
of the effectiveness of AMR on the quality of MHD turbu- 
lence simulations. By varying the refinement criteria, the 
base grid resolution, and the number of refinement levels, 
we find that the quality of the turbulence statistics — in 
particular, the spectral index of the velocity power spec- 
trum and the extent of the inertial range — is more closely 
related to the average refinement coverage than to the 
maximum level of refinement. Analysis of the density 
PDFs shows that, with our refinement criteria, AMR is 
particularly powerful for simulations in which interest is 
focused on the regions of highest density: it does not 
capture the regions of the lowest density as well as a 
unigrid simulation, but it does capture th e PDF of the 
high-density regions quite accurately (e.g. iCollins et al.l 
[201 111 . Our result for the dissipation coefficient for tur- 
bulence with = 10 and /3o = 0.1 is e = 0.48 ± 0.01. 
This is consistent with other numerical studies of MHD 
turbulence with unigrid codes, which find that e ~ 0.5. 
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Fig. 1. — One of the iRvu &: Jonej II1995I) shock-tube test results using the ORION2 MHD module. The top row, from left to right, shows 
the density (p), pressure (P), total energy density (-Etot) along the shock tube, which is resolved by 512 cells. The second row shows the 
three components of velocity and the bottom row shows the three components of magnetic field. In these panels, from left to right, we can 
see the fast shock, slow rarefaction (apparent in the plots of Vy and By), contact discontinuity, slow shock, and fast shock. 
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Fig. 2. — IBrio fc Wul II1988I ) shock tube test results. The top row, from left to right, shows the density (p), pressure (P), and velocity 
along the shock tube (fx). The second row shows the y-component of the velocity (vy), the temperature (T), and the y-component of the 
magnetic field {By). The shock tube is resolved by 800 cells. A compound wave composed of an Alfven and a slow wave is seen as the 
spike near the middle of the figure. 
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(c) ^ 




Fig. 3. — Vertical magnetic field-loop advection test on a unigrid with resolution A'^ = 64 (a) at the beginning of the advection and 
(b) after 2 cycles of advection. The inclined magnetic field loop, which is rotated by 63.4° about the y-axis from the vertical, (c) at the 
beginning of the advection and (d) after 2 cycles of advection. Both field loops are advected with velocity (2,1,1). All four panels are at 
z=0.5. The field loops maintain their shapes accurately after 2 cycles of advection. (A color version of this figure is available in the online 
journal.) 
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Fig. 4. — Time evolution of (a) the mean-squared magnetic field, (B^), normalized by the initial value at time = and (b) the normalized 
error, (|B3|)/_Bo, for the inclined magnetic field-loop advection tests. The three unigrid tests are shown as thin dot-dashed (N = 32), thin 
dashed (A'' = 64), and thin solid (A'' = 128) curves. The dotted curve (which nearly overlaps the N = 128 unigrid test result) is from the 
advection test with AMR at A'^ = 64 resolution on the base level and 1 level of refinement (see Figure [6]l. The thick dashed curve is from 
the advection test with Fixed Mesh Refinement (FMR) at A'^ = 64 resolution at the base grid and 1 level of refine ment on a fixed fine grid 
smaller than the base grid (see FigureO. Both the AMR and the FMR tests use a refinement ratio n^cf = 2. See i]3.2l for discussion of the 
tests. 
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Fig. 5. — Volume rendering of the inclined field loop from the 3D N = 64 unigrid model after two cycles of advection along the box 
diagonal. (A color version of this figure is available in the online journal.) 
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Fig. 6. — Inclined magnetic field-loop advection test with AMR at the beginning of the advection (left panel) and after 2 cycles of 
advection (right panel). The blocks surrounding the loop show the location of the fine grids. Although the base grid has a resolution of 
N = 64, the loop is always refined at a resolution of N = 128. The time evolution of the normalized (-B^) and (I-B3I) is almost the same as 
that of a unigrid A'^ = 128 test (Fig. IJ}. (A color version of this figure is available in the online journal.) 
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M (d) 




Fig. 7. — Inclined magnetic field-loop advection test with fixed mesh refinement at (a) time = 0, (b) time = 0.25, (c) time = 0.5, and (d) 
time = 2.0. The base grid resolution is A' = 64 and the fine grid is equivalent to = 128. When the field loop crosses the coarse-fine grid 
boundary, the loop diffuses to the resolution of = 64 but maintains its shape accurately after 2 cycles of advection. See Figure|4]for the 
time evolution of the normalized {B^) and (I-B3I) for this fixed mesh test. (A color version of this figure is available in the online journal.) 




k 

Fig. 8. — Compensated velocity power spectra of models 1, 3, 4, 6, and 10. The spectra are compensated by fc^'*^, the spectral index 
of the unig rid 512^ Model 10. Models 3, 4, and 6 are simulations with AMR. Models 1 and 10 have just the base grid at 128"^ and 512'', 
respectively. The sp ectr a are labeled by their model numbers. The horizontal line is provided to aid comparison of the inertial ranges of 
all the models. See i|4.3l for a discussion of how refinement coverage affects the inertial range. (A color version of this figure is available in 
the online journal.) 
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Fig. 9. — Magnetic field power spectra of models 1, 3, 4, 6, and 10 (with model numbers shown in the same order as the curves, starting 
from bottom) without compensation. The spectra show a similar convergence behavior with resolution as the velocity spectra in Figure [81 
(A color version of this figure is available in the online journal.) 




Fig. 10. — Time-averaged velocity power spectra of Model 10. The solenoidal Pv^ik) (dot-dashed) and the dilatational P„^(fc) (dashed) 
components of the velocity spectrum are also shown. The time-averaged fraction of dilatation component Xc{k) = 0.26 ± 0.02. The 
solenoidal component dominates as a result of purely solenoidal driving of the turbulence. 
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log p log p 

Fig. 11. — Comparison of density PDFs. Panel (a) shows unigrid Model 1 (dot-dash red line), AMR Model 4 (dashed blue line), and 
unigrid Model 10 (solid black line). The AMR model matches the high-resolution unigrid Model 10 very well at high densities, but it is 
similar to the low resolution unigrid Model 1 at low densities. Panel (b) shows AMR Model 2 (dot-dash red line), AMR Model 3 (dashed 
blue line), and unigrid Model 10 (solid black line). Model 2 has a very low volume coverage at refinement level 2 and as a result the high 
density end deviates more from the high-resolution unigrid model than does Model 3. (A color version of this figure is available in the 
online journal.) 
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Fig. 12. — Turbulent dissipation rate coefficient, e, of tfie 10 MHD turbulence models listed in Table 1. For (R) > 350, e agrees with the 
value in the 512^ unigrid simulation, 0.476 ± 0.011. Model 9 has (R) ~ 624 and e = 0.467 ± 0.013. Model 11 (solid circle) has a somewhat 
smaller value of e than Model 3, which has the same value of the index of the velocity power spectrum. 
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Fig. 13. — Normalized enhancement of the magnetic field strength, (|A_B|)/So, of the 11 MHD turbulence models listed in Table 1. 
Model 11 {Ma,o = 1; solid circle) has a smaller enhancement in the magnetic field strength than the 10 models with Ma,o = V5. 



TABLE 1 

Velocity Power Spectral Indexes, Inertial Ranges, and Refinement 



Model 


base 


Refinement 


Refinement 


Shear flow 


Spectral 


Inertial 


Refinement 


(R) 


Normalized 




grid 


levels 


threshold 


refinement* 


index (n) 


range'' 


coverage 


= (%) 


CPU time'^ 
















e = 1 


e = 2 






1 


128 









1.75 ±0.06 


4 r. 


- 13 






128 


3.81(-3) 


2 


128 


2 


3.25 


no 


1.76 ±0.07 


4 - 


- 13 


6 


1.3 


139 


5.17(-2) 


3 


128 


2 


2.75 


no 


1.65 ± 0.06 


4 - 


- 14 


36 


9 


197 


4.87(-l) 


4 


128 


2 


2.5 


no 


1.61 ± 0.03 


4 - 


- 17 


65 


17 


255 


6.30(-l) 


5 


128 


2 


2.5 


yes 


1.58 ±0.03 


4 - 


- 17 


70 


18 


264 


6.58(-l) 


6 


128 


1" 


2.5 


yes 


1.48 ±0.04 


4 r. 


^ 21 


76 




420 


2.46 


7 


256 


1 


2.5 


no 


1.46 ± 0.04 


4 - 


- 18 


19 




304 


2.41 


8 


256 


2 


2.5 


yes 


1.48 ± 0.04 


4 r. 


- 21 


31 


6.5 


369 


4.71 


9 


256 


2 


2.3 


yes 


1.46 ±0.03 


4 - 


- 25 


58 


17 


625 


9.47 


10 


512 









1.42 ± 0.02 


4 - 


- 26 






512 


1.0 


llf 


128 


2 


2.75 


yes 


1.65 ±0.06 


4 - 


- 14 


40 


10.4 


206 


4.05(-l) 



^ When the shear flow refinement eriterion is not included, refinement is determined by only the total pressure jump. 
^ The inertial range extends from k — 4 to k — fcmax- 

i — i stands for volume coverage at level i. The standard deviation of fluctuations in refinement coverage is ^ 3% 
^ The total CPU time is normalized by the total CPU time of the Model 10 unigrid run. 

Refinement ratio n^^f — 4, corresponding to a maximum resolution of 512'^. 
^ Model 11 is a stress test of the code and has A^rms 17.32, /3o = 0.00667 and A4a,o — li in contrast to models 1-10, which have Alrms = 10, 
/3o = 0.1 and A^a,o = "s/S- 



